setwd("/Volumes/HD_2/Documents/lab/project/OIS/450Karray/1117_Amir/PROCESSED_DATA_FILES")
####test
amir_beta<-read.table("1117_AmirEden_CORRECTED_BETA_VALUES.txt",sep="\t",header=T,skip=22)
beta_test<-data.frame(pBabe1=amir_beta[,10],pBabe2=amir_beta[,11],pBabe3=amir_beta[,12],Ras1=amir_beta[,13],Ras2=amir_beta[,14],Ras3=amir_beta[,15])

p <- ggplot(beta_test, aes(pBabe1, pBabe2))
p + geom_point(alpha=1/10,size=1.2)

p <- ggplot(beta_test, aes(pBabe1, Ras1))
p + geom_point(alpha=1/100,size=1)

p <- ggplot(beta_test, aes(Ras1, Ras2))
p + geom_point(alpha=1/50,size=1.2)

p <- ggplot(beta_test, aes(pBabe1, pBabe3))
p + geom_point(alpha=1/50,size=1.2)

p <- ggplot(beta_test, aes(pBabe2, pBabe3))
p + geom_point(alpha=1/50,size=1.2)

p <- ggplot(beta_test, aes(Ras1, Ras3))
p + geom_point(alpha=1/100,size=1)

p <- ggplot(beta_test, aes(pBabe2, Ras2))
p + geom_point(alpha=0.5,size=1)

p <- ggplot(beta_test, aes(pBabe2, Ras1))
p + geom_point(alpha=1/50,size=1.2)

####Beginning!!

### load data
mainfest<-read.delim("/Volumes/HD_2/Documents/lab/project/OIS/450Karray/GPL13534-10305.txt",sep="\t",header=T,skip=37)
amir_beta<-read.table("1117_AmirEden_CORRECTED_BETA_VALUES.txt",sep="\t",header=T,skip=22)
mask_probes<-read.table("/Volumes/HD_2/Documents/lab/project/OIS/450Karray/MaskedProbes.450k.txt",sep="\t",header=F)
mainfest.sort<-mainfest[order(mainfest[,1]),]


detection_p<-read.table("1117_AmirEden_P-VALUES.txt",sep="\t",header=T,skip=22)
detection_p<-detection_p[,10:15]
beta.manifest<-cbind(mainfest.sort,amir_beta[,10:15],detection_p)

### mask probes with SNP, repeats, any NA value among 6 samples..
beta.manifest.masked<-beta.manifest[!(beta.manifest$ID %in% mask_probes[,1]),]
beta.manifest.masked.noNa<-beta.manifest.masked[rowSums(is.na(beta.manifest.masked[,38:43]))==0,]
colnames(beta.manifest.masked.noNa)[38:43]<-c("pBabe1","pBabe2","pBabe3","Ras1","Ras2","Ras3")

beta.manifest.masked.noNa.CpgProbes<-beta.manifest.masked.noNa[grep("cg", beta.manifest.masked.noNa$ID),]
beta.manifest.masked.noNa.CphProbes<-beta.manifest.masked.noNa[grep("ch", beta.manifest.masked.noNa$ID),]

samples.OIS<-c("pBabe1","pBabe2","pBabe3","Ras1","Ras2","Ras3")

## detection p value shows that pBabe1 sample is an outlier.
dim(beta.manifest.masked.noNa[rowSums(beta.manifest.masked.noNa[,44:49] <0.05)==0,])
dim(beta.manifest.masked[beta.manifest.masked[,44] <0.05,])


### scatter plot
##multiple pdf generated, only OIS samples
beta.manifest.masked.noNa.OIS<-data.frame(beta.manifest.masked.noNa[,38:43],row.names=rownames(beta.manifest.masked.noNa[,1]))
for(i in c(1:dim(beta.manifest.masked.noNa.OIS)[2])){
	for(j in c(1:dim(beta.manifest.masked.noNa.OIS)[2])){
		if(i<j){
			#filename=paste(sub(".AVG_BETA", "", colnames(beta.manifest.masked.noNa.OIS)[i]),"_vs_",sub(".AVG_BETA", "", colnames(beta.manifest.masked.noNa.OIS)[j]),".jpeg",sep="")
			filename=paste(samples.OIS[i],"_vs_", samples.OIS[j],".2by2.AVG_BETA.jpeg",sep="")
			jpeg(filename)
			p <- ggplot(beta.manifest.masked.noNa.OIS, aes(beta.manifest.masked.noNa.OIS[,i], beta.manifest.masked.noNa.OIS[,j])) + geom_point(alpha=1/100,size=1) + labs(list(title=filename, x = samples.OIS[i], y = samples.OIS[j]))
			print(p)
			dev.off()
		}
		
	}
}

##multiple pdf generated, only OIS samples, in smoothScatter
beta_test
smoothScatter(beta_test$pBabe2, beta_test$Ras1, nbin = 1000, xlab = "pBabe2", ylab = "Ras1", xlim=c(0,1), ylim=c(0,1))

beta.manifest.masked.noNa.OIS<-data.frame(beta.manifest.masked.noNa[,38:43],row.names=rownames(beta.manifest.masked.noNa[,1]))
for(i in c(1:dim(beta.manifest.masked.noNa.OIS)[2])){
	for(j in c(1:dim(beta.manifest.masked.noNa.OIS)[2])){
		if(i<j){
			#filename=paste(sub(".AVG_BETA", "", colnames(beta.manifest.masked.noNa.OIS)[i]),"_vs_",sub(".AVG_BETA", "", colnames(beta.manifest.masked.noNa.OIS)[j]),".jpeg",sep="")
			filename=paste(samples.OIS[i],"_vs_", samples.OIS[j],".2by2.AVG_BETA.smoothScatter.jpeg",sep="")
			title=paste(samples.OIS[i],"_vs_", samples.OIS[j],sep="")
			jpeg(filename)
			smoothScatter(beta.manifest.masked.noNa.OIS[,i], beta.manifest.masked.noNa.OIS[,j], nbin = 1000, xlab = samples.OIS[i], ylab = samples.OIS[j], xlim=c(0,1), ylim=c(0,1), font.lab=2, font=2, main=title)
			abline(0,1)
			dev.off()
		}
		
	}
}

##multiple pdf generated, only OIS samples, alpha=1 for outlier
beta.manifest.masked.noNa.OIS<-data.frame(beta.manifest.masked.noNa[,10:15],row.names=rownames(beta.manifest.masked.noNa[,1]))
for(i in c(1:dim(beta.manifest.masked.noNa.OIS)[2])){
	for(j in c(1:dim(beta.manifest.masked.noNa.OIS)[2])){
		if(i<j){
			#filename=paste(sub(".AVG_BETA", "", colnames(beta.manifest.masked.noNa.OIS)[i]),"_vs_",sub(".AVG_BETA", "", colnames(beta.manifest.masked.noNa.OIS)[j]),".outlier.jpeg",sep="")
			filename=paste(samples.OIS[i],"_vs_", samples.OIS[j],".2by2.AVG_BETA.outlier.jpeg",sep="")
			jpeg(filename)
			p <- ggplot(beta.manifest.masked.noNa.OIS, aes(beta.manifest.masked.noNa.OIS[,i], beta.manifest.masked.noNa.OIS[,j])) + geom_point(alpha=1,size=1) + labs(list(title=filename, x = samples.OIS[i], y = samples.OIS[j]))
			print(p)
			dev.off()
		}
		
	}
}


##corelation
cor_matrix=""
for(i in c(1:dim(beta.manifest.masked.noNa.OIS)[2])){
	for(j in c(1:dim(beta.manifest.masked.noNa.OIS)[2])){
		if(i<j){
			cor_matrix=c(cor_matrix,cor(beta.manifest.masked.noNa.OIS[,i], beta.manifest.masked.noNa.OIS[,j],use= "complete.obs", method = "pearson"))
			
		}
		
	}
}

### select most variant probes (range, variance to do histogram firstly), do distance matrix 
variance<-apply(beta.manifest.masked.noNa[,38:43],1,var)
ran<-function(x){return(max(x)-min(x))}
range<-apply(beta.manifest.masked.noNa[,38:43],1,ran)

summary(variance)
summary(range)
hist(variance,xlim=c(0,0.03),breaks=seq(0,0.13,by=0.0001))
hist(range,xlim=c(0,1),breaks=seq(0,1,by=0.01))

beta.manifest.masked.noNa.varRange<-cbind(beta.manifest.masked.noNa,variance,range)
#beta.manifest.masked.noNa.varSig<-beta.manifest.masked.noNa.varRange[beta.manifest.masked.noNa.varRange[,44]>0.0005,]
#beta.manifest.masked.noNa.varSig0.01<-beta.manifest.masked.noNa.varRange[beta.manifest.masked.noNa.varRange[,44]>0.014,] 3242
beta.manifest.masked.noNa.varSigTop3000<-beta.manifest.masked.noNa.varRange[order(beta.manifest.masked.noNa.varRange[,44],decreasing=T),]
beta.manifest.masked.noNa.varSigTop3000<-beta.manifest.masked.noNa.varSigTop3000[1:3000,]
beta.manifest.masked.noNa.rangeSig<-beta.manifest.masked.noNa.varRange[beta.manifest.masked.noNa.varRange[,45]>0.1,]
#beta.manifest.masked.noNa.rangeSig0.3<-beta.manifest.masked.noNa.varRange[beta.manifest.masked.noNa.varRange[,45]>0.3,]

cl<-hclust(dist(t(beta.manifest.masked.noNa.varSigTop3000[,38:43])), method = "average", members=NULL)
plot(cl,main="Distance between Samples\n(Top 3000 most variant probes)",ylab="Distance",xlab="Sample names")

cl<-hclust(dist(t(beta.manifest.masked.noNa.varRange[beta.manifest.masked.noNa.varRange[,44]>0.0005,][,38:43])), method = "average", members=NULL)
plot(cl,main="Distance between Samples\n(Variance>0.0005, 196052 probes)",ylab="Distance",xlab="Sample names")

cl<-hclust(dist(t(beta.manifest.masked.noNa.varRange[beta.manifest.masked.noNa.varRange[,45]>0.05,][,38:43])), method = "average", members=NULL)
plot(cl,main="Distance between Samples\n(Range>0.05, 214882 probes)",ylab="Distance",xlab="Sample names")

### t-test , calculate q value


pvalue.varSigTop3000=-1
for(i in c(1:length(beta.manifest.masked.noNa.varSigTop3000[,1]))){
	p<-t.test(as.numeric(beta.manifest.masked.noNa.varSigTop3000[i,38:40]), as.numeric(beta.manifest.masked.noNa.varSigTop3000[i,41:43]), alternative = "two.sided",exact=T, paired = F, na.rm=T)
	pvalue.varSigTop3000<-c(pvalue.varSigTop3000,as.numeric(p$p.value))
}
pvalue.varSigTop3000<-pvalue.varSigTop3000[2:length(pvalue.varSigTop3000)]


library(qvalue)
qvalue.varSig<-qvalue(pvalue.varSigTop3000)

qsummary(qvalue.varSig)



### limmar calculate adjusted p value (cant use it..., varSig0.0005 is ok.. there are 77 p<0.05 probes)
design <- cbind(Grp1=1,Grp2vs1=c(0,0,0,1,1,1))
fit <- lmFit(beta.manifest.masked.noNa.varSig[,38:43],design)
fit <- eBayes(fit)
topTable(fit,coef=2,adjust.method="BH")
summary(topTable(fit,coef=2,adjust.method="BH",number=1000000)$adj.P.Val)
Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.08476 0.38830 0.65730 0.61880 0.85950 1.00000





### do vocano plot, select most possible changed probes
betaDiff<-rowMeans(beta.manifest.masked.noNa.varSigTop3000[,41:43],na.rm=T)-rowMeans(beta.manifest.masked.noNa.varSigTop3000[,38:40],na.rm=T)
beta.manifest.masked.noNa.varSigTop3000<-cbind(beta.manifest.masked.noNa.varSigTop3000,qvalue.varSig$qvalues,betaDiff)

dim(beta.manifest.masked.noNa.varSigTop3000[beta.manifest.masked.noNa.varSigTop3000[,47]< -0.2|beta.manifest.masked.noNa.varSigTop3000[,47]>0.2,])
dim(beta.manifest.masked.noNa.varSigTop3000[beta.manifest.masked.noNa.varSigTop3000[,47]< -0.2,])
dim(beta.manifest.masked.noNa.varSigTop3000[beta.manifest.masked.noNa.varSigTop3000[,47]>0.2,])

beta.manifest.masked.noNa.varSigTop3000.sigPandBeta<-beta.manifest.masked.noNa.varSigTop3000[beta.manifest.masked.noNa.varSigTop3000[,46]< 0.05 &(beta.manifest.masked.noNa.varSigTop3000[,47]< -0.2|beta.manifest.masked.noNa.varSigTop3000[,47]>0.2),]
beta.manifest.masked.noNa.varSigTop3000.sigPandBeta.hyper<-beta.manifest.masked.noNa.varSigTop3000[beta.manifest.masked.noNa.varSigTop3000[,46]< 0.05 & beta.manifest.masked.noNa.varSigTop3000[,47]>0.2,]
beta.manifest.masked.noNa.varSigTop3000.sigPandBeta.hypo<-beta.manifest.masked.noNa.varSigTop3000[beta.manifest.masked.noNa.varSigTop3000[,46]< 0.05 & beta.manifest.masked.noNa.varSigTop3000[,47]< -0.2,]


vocanoData<-data.frame(fdr=-log10(beta.manifest.masked.noNa.varSigTop3000[,46]),betaDiff=beta.manifest.masked.noNa.varSigTop3000[,47])
jpeg("Vocano_plot.MeanBeta_vs_qValue.jpg")

p <- ggplot(vocanoData, aes(vocanoData$betaDiff, vocanoData$fdr)) 
p + geom_point(alpha=1,size=1.5) + labs(list(x = "Mean Beta differences", y = "-log10 Adjusted P value")) + geom_hline(yintercept = -log10(0.05)) + geom_vline(xintercept = -0.2) + geom_vline(xintercept = 0.2)
dev.off()

### use these probes, do genomic region enrichment (genes, CGI, Enhancer)
#CGI

dim(beta.manifest.masked.noNa.varSigTop3000.sigPandBeta[grep("Island", beta.manifest.masked.noNa.varSigTop3000.sigPandBeta$Relation_to_UCSC_CpG_Island),])

phyper(2,142173,394029-142173,46,lower.tail = F)


### use these probes, do clustering and heatmap
library(colorRamps)
library(gplots)
mat<-cbind(as.numeric(unlist(beta.manifest.masked.noNa.varSigTop3000.sigPandBeta[,38])),as.numeric(unlist(beta.manifest.masked.noNa.varSigTop3000.sigPandBeta[,39])),as.numeric(unlist(beta.manifest.masked.noNa.varSigTop3000.sigPandBeta[,40])),as.numeric(unlist(beta.manifest.masked.noNa.varSigTop3000.sigPandBeta[,41])),as.numeric(unlist(beta.manifest.masked.noNa.varSigTop3000.sigPandBeta[,42])),as.numeric(unlist(beta.manifest.masked.noNa.varSigTop3000.sigPandBeta[,43])))
pdf("heatmap.beta.OIS.pdf")
heatmap.2(mat,
		Rowv=TRUE,
		Colv=TRUE,
		dendrogram= c("both"),
		distfun = dist,
		hclustfun = hclust,
		na.rm=TRUE,
#xlab = "Sample", ylab = "ASM blocks",main = "ASM Blocks -log10 P value",
		labCol = c("pBabe1","pBabe2","pBabe3","Ras1","Ras2","Ras3"),
		labRow = NULL,
#angle=45,
		key=TRUE,
		keysize=1,
		trace="none",
		density.info=c("none"),
		margins=c(10, 8),
#breaks=c(-10:10),
#symbreaks=TRUE,
		na.color=par("bg"),
#col=brewer.pal(10,"PiYG")
		#col=redgreen(75),
		col=blue2yellow(75),
)
dev.off()

##all of top 3000 most variant probes
mat<-cbind(as.numeric(unlist(beta.manifest.masked.noNa.varSigTop3000[,38])),as.numeric(unlist(beta.manifest.masked.noNa.varSigTop3000[,39])),as.numeric(unlist(beta.manifest.masked.noNa.varSigTop3000[,40])),as.numeric(unlist(beta.manifest.masked.noNa.varSigTop3000[,41])),as.numeric(unlist(beta.manifest.masked.noNa.varSigTop3000[,42])),as.numeric(unlist(beta.manifest.masked.noNa.varSigTop3000[,43])))
pdf("heatmap.beta.OIS.allMostVaraintProbes.pdf")
heatmap.2(mat,
		Rowv=TRUE,
		Colv=TRUE,
		dendrogram= c("both"),
		distfun = dist,
		hclustfun = hclust,
		na.rm=TRUE,
#xlab = "Sample", ylab = "ASM blocks",main = "ASM Blocks -log10 P value",
		labCol = c("pBabe1","pBabe2","pBabe3","Ras1","Ras2","Ras3"),
		labRow = NULL,
#angle=45,
		key=TRUE,
		keysize=1,
		trace="none",
		density.info=c("none"),
		margins=c(10, 8),
#breaks=c(-10:10),
#symbreaks=TRUE,
		na.color=par("bg"),
#col=brewer.pal(10,"PiYG")
		#col=redgreen(75),
		col=blue2yellow(75),
)
dev.off()
